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We show that any non-relativistic quantum A'^-body dynamics problem with pairwise interac- 
tions can be exactly reformulated in terms of A'^ well behaved 1-body stochastic density equations. 
Specifically, the time evolving A''-body density matrix is written as an average of tensor products 
of 1-body densities each of which obeys a stochastic evolution equation. Such decompositions can 
be constructed for any mixture of fermions and bosons. The evolution equations for the 1-body 
densities preserve norm, Hermiticity and positivity. 
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While the dynamical laws underlying much of chem- 
istry and physics are known, these equations can only 
rarely be solved analytically, and in fact cannot even be 
solved numerically except when the number of partici- 
pating particles is very small. This is a major difficulty 
for theoretical chemistry and solid state physics where 
systems of many electrons must be solved in order to 
determine electronic structure. Approximation methods 
such as density functional theory and Hartree-Fock and 
its generalizations allow estimates of ground state ener- 
gies and other properties, but limited accuracy and gen- 
erality is inherent in such approaches. Other more for- 
mal computational approaches exist, which are based on 
Lanczos type algorithms, configuration interaction^ and 
coupled-cluster [31 schemes, which converge in principle 
to exact solutions in appropriate limits. The difficulty of 
such methods stems from the practical necessity of rep- 
resenting the wavefunction in a basis of 1-body states[3j. 
That is, if m basis functions are required for each body 
then complex numbers will be needed to represent 
the wavefunction. This exponential growth of computa- 
tional complexity with N is characteristic of almost all 
numerical strategies which converge to exact solutions as 
TO — !■ oo. Very recently several classes of computational 
methods have been discovered which appear to circum- 
vent this prohibitive scalingjj 0,0,0, lA]- 

The methods of interest here are those based on wave 
equations and stochastic representation of the pairwise 
interaction 0,0,0|- Polynomially scaling schemes for evo- 
lution of systems of identical bosons and fermions 
have been developed which decompose the state vector 
as an average of solutions of 1-body stochastic wave equa- 
tions. These methods avoid explicit construction or stor- 
age of the state vector and require storage which scales 
only linearly with N . Actual computational costs scale 
as some small power of N . Unphysical properties of such 
deconipositions, such as exponential growth of 1-body 
norm0,0|, have been eliminated by rederivation of the 
evolution equations'^ from a variational principle|3|. In 
addition, considerable progress has been made in devel- 
oping numerical strategies for solving the 1-body stochas- 
tic evolution equations!^. These developments suggest 



that stochastic reformulations of the iV-body problem 
may ultimately lead to more effective computational 
methods for quantum chemistry and physics. 

An interesting open question is whether it is possi- 
ble to reformulate all non-relativistic many body prob- 
lems using such stochastic decompositions? Existing de- 
compositions work only for systems of identical fermions 
or systems of identical bosons. Mixtures are not al- 
lowed. These restrictions on composition seem artificial 
and limit application. Fundamental physical systems as 
simple as heteronuclear molecules cannot be evolved with 
these methods. In this manuscript we develop a general 
unified theory applicable to any non-relativistic A^-body 
system of pairwise interacting distinguishable and indis- 
tinguishable quantum particles. This is achieved by ex- 
pressing the A'^-body density matrix as an average of ten- 
sor products of A^ 1-body densities each of which obeys 
a stochastic evolution equation. We also show that the 
stochastic evolution equations are surprisingly well be- 
haved with preservation of 1-body norm, Hermiticity and 
positivity. Our variational derivation also shows that the 
decomposition with these favorable properties is unique. 
Finally, we provide a simple recipe for recovery of the 
state vector from the A^-body density matrix. 

Consider then the general Hamiltonian for a system of 
particles with pairwise interactions 

N N-1 N 

nN = Y.Hk + Y. E (1) 

fc=l fc=l l=k+l 

where Hk denotes the single body Hamiltonian for par- 
ticle k and Vk^i is a general interaction potential for par- 
ticles k and I. Without loss of generality we write the 
general 2-body interaction Vk^i as a sum of products of 
1-body operators 

Vk,i^j2''sOlOt (2) 

where p may be arbitrarily large. [Note that for a 
Coulomb interaction p is in principle infinite, but in prac- 
tice restriction to a finite basis of 1-body states leads to 
finite p.] This decomposition can always be chosen such 
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that the dimensionless 1-body operators are Hermi- 
tian. The real coefficients Ug have units of energy and 
may be positive or negative. 

Leaving issues of symmetry due to indistinguishability 
for later, consider an initial iV-body density matrix which 
is a product of 1-body densities. Exact dynamics of the 
N-body density p^{t) will be governed by the Liouville- 
von Neumann equation. We will look for a representation 
of the time evolved density as an average tensor product 
of 1-body densities. Specifically, we will derive a decom- 
position of the form 



(3) 



where Pkit) for k = 1, . . . , N are intended to be one-body 
density operators which obey stochastic evolution equa- 
tions yet to be determined, and M[. . .] denotes an average 
over different realizations of the stochastic processes. Our 
decision to decompose the density matrix rather than the 
state vector is motivated by a proof 10] that no such gen- 
eral decomposition conserves 1-body norm. 

Introducing a set of complex Wiener processes af. ; (t) 
the most general form of the stochastic evolution equa- 
tions is 

p N 
dpk{t) - Vk{t)dt + Y,ul(t)Y,dali(t) 



1=1 

N 



s=l l^k 



(4) 



where Vk{t) and u%{t) are unknown operators to be de- 
termined, and where Vk{t) must be Hermitian in order 
that Hermiticity be preserved. The complex stochastic 
differentials are of the form (/i -I- iv)^fdi where 

p,, V are normally distributed with zero mean and unit 
variance. Additionally, we require that Trfc{ufe} = and 
Trfe{M|} = so that trace norm will be conserved. If we 
employ Ito calculus0| then 
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dpkUpi+Y. J2 '^PkdPi n '°'"] 



fc=l 



l^k 



k=l l=k+l 



ra^k.l 



and we will require that this equal —i[Ti.N,p ]dt so that 
our decomposition is exact (note that in our units h = 1). 
That is, we require that the average tensor product equal 
the exact A'^-body density matrix p^ which obeys the 
Liouville-von Neumann equation dp^ = —iCp^dt, L be- 
ing the Liouville super-operator L, = [TYat, . ]. Substitu- 
tion of 0) into © reveals that we must have terms in 
which are proportional to in order that the pair- 
wise coupling terms in Tijv be reproduced. In addition 
to avoid unwanted terms it will be necessary to require 



that da 



Lk 



daf. *i and that the Wiener processes be oth- 



erwise independent. Standard rules of Ito calculus|llj 
then determine the following further relations 



dal*idal, 



5k,k'^i,i' 



(6) 



among stochastic differentials. 

To find explicit forms for Vk and by trial and error 
would be difficult. The task is greatly simplified by use 
of a variational principle^. Defining a Hilbert-Schmidt 
norm ||a:|P = {x\x) and inner product {x\y) — Tr{a;^y} 



(see 113) we may express the mean square error in a 
decomposition as a functional 

T ^\dp^ + iCp^dt\\^ (7) 

which can be minimized by variation of and u^.. Al- 
lowing all independent variations 6vk, Su^, and Su^i^\ 
substituting Eqs. ^ and (^J into Eq. (0, keeping only 
terms of order dt and dropping dt as a factor, then gives 
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= ST = TT{{Y^5vkY[ 
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fc=l l=k+l s=l 
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N-1 N 



Vk 



'Up^ + E E Ek-^ + 



fc=l l^k k=l l=k+l s=l 

i Hnp"" - I p'^Un)} 



N 

n ^™ 

(8) 



where (5['u^uf -I- -u^^w^^] = Sulu^ + ulSuf + Sul'^u'^'^ + 
u^f.^Sui^ . The solution of Eq. © still involves tedious 
algebra and hence the complete derivation will be given 
elsewhere [ifij. One finds that 

p N 

dpkit) = -i[Hk,Pk{t)]dt~iy^y^L0sO"{t)[oiPk{t)]dt 



s=l l^k 



p 



N 



+ >^V^(0^-0,(i)K(t) yjalAt) 



s=l 



l^k 



N 



+ T(v^srpkitm-oiit))ydaim (9) 



s=l 



l^k 



where Oj,(t) = Trk{Olpk{t)}. A detailed examina- 
tion of this derivation indicates that these equations are 
uniaue[lo|. We now focus on the properties of these evo- 
lution equations for the stochastic one-body density op- 
erators. 

Proof of exactness: Substituting Eq. jnj into us- 
ing Eqs. ©I recognising that the differentials are statis- 
tically independent of the densities at time t, and using 
the fact that Af [c?a|. J = we have 
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Hk,Pk]ll 



Pi 



l^k 



N N p N 
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(10) 
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N-l N 



N 



k=l l=k+l s=l m-^k.l 
N-l N p N 



where we have used the fact that Tr{^i3} = Tv{BA} 
and the definition of O^. 

Proof of positivity: Surprisingly, the one-body density 
operators remain positive semi-definite. To see this recall 



that the Hermitian property of the stochastic one-body 



l=k+l s=l 



By definition p^ = ^[Pk]\f-Lk Pi\ ^ M[pfep/ Hm^^fc,; Pm]> 
and so after some algebra it follows that 



,N 



density operators permits a spectral decomposition of the 
form 



N 



(13) 



dp 



N 



-iJ2[Hk,p'']dt 



k=l 
N-l N p 



(11) 



- *E E J2''^[oiot,p'']dt 

fc=l l=k+l s=l 

N N p 

- *EEE^^^;[o^,p^]dt 

fe=l l^k s=l 
N-l N p 

+ *E E E^4o:[o^,p^]+o:,[of,p^]} dt. 

k=l l=k+l s=l 

Finally, using the following relation of the summation 
operators in Eq. lfTT|) 



where £l{t) = \k, j; (t)) {k, j; {t)\ are the spectral pro- 
jection operators for the natural orbitals and u'^(i) are 
their occupation probabilities (i.e. pk{t)\k, j; {t)) — 
'wl,{t)\k,j;{t))). The question is whether the u;^(i) are 
positive for all time. Hence, consider the following 
lemma: 

Lemma: The decomposition H13|l satisfying Eqs. 
implies that 

p N 
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l^k 
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EE[-]^[ 

k=l l^k 



N-l N 

]. = E E [-i^t-]' 

k=l l=k+l 



Y,iV^.rm),~oi)wiit)j2dai} (u) 

s=l l^k 

p N 



dSiit) 



Hk.Slit) 



dt - i^^uJsO^i 

s=l l^k 



dt 



the last two terms in Eq. cancel out and we obtain 
the Liouville-von Neumann equation 
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N-l N p 



E^^ + E E E-^of-of, 



L fc=l k=l l=k+l s=l 

N ■ 



^ -i[nN,p'^] dt 



dt 
(12) 



- (N- l)j2^sm), - Ol){Ol - {Oi),)Si{t) dt 

s=l 

- iN-i)j2^s£i{t)m),-oi)ioi - {01),) dt 

p N 

+ J2 V^siOl - {Oi)i)£i{t) dali 



which proves that the decomposition is exact. 

Proof of Hermiticity: To show pk (i) = pj, (i) take the 
adjoint of Eq. 0. Since this results in the same equation 
with the same initial condition it follows that pkit) is 
Hermitian. 

Proof of norm conservation: In order to prove that 
trace-norm is preserved we trace over Eq. (O which 
yields 

Ti'kidpk} = -iTTk{[Hk, Pk]} dt 

p N 

- «EE^^^;Tr,{[0^.,pfc]}dt 

s=l l^k 

p N 

+ J2V^sTTk{{Ol~Ol)pk}J2dali 



s=l 

p 



l^k 



N 



j2{V^s)*simoi - {oi),) J2 rf« 

s=l l^k 



s * 
k,l- 



(15) 



Here (O^), ^ Tr,{0^£^(t)}. 
Proof: From H13|l it follows that 

oo oo 

dpk = J2dwi{t)£l{t) + Y,wi{t)d£l{t) 



OO 



+ Y.dwiit)d£lit) 



(16) 



l^k 



N 



and substituting Eqs. (|14|) and (|15|) and using the rela- 
tions ijHJl then gives 

Finally, from Eq. (|14|l it readily follows that the occu- 
pation probabilities must satisfy the relation 



Y,{V^s)*TTk{pkiOt, - Ol)} J2 da 

s=l l^k 





s * 

k,l 



p N 

s=l l^k 
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for T sufBciently large. 



/ m),{t')-oUt')rdt'} 



> 



(17) 



which shows that the occupation probabihties are indeed 
positive. It follows that the 1-body densities are positive 
semi-definite. 

Finally, we outline a simple method to extract the wave 
vector from the density operator. Suppose the initial 
state for the iV-body system is a pure state vlf^(O). Now 
choose a set of arbitrary states \ik), one for each particle. 
Then construct 

= M [piit)\h) ® P2m2) ® • • • ® PN{t)\lN)]{lS) 

and normalize to obtain |l>^(t))/||'l^(i)||. 
This vector now differs from the exact wave function by 
a possible phase factor: that is 



= |*^^(t))exp[ze(t)] 



(19) 



with d|*^(t)) = -iHN\'^^{t))dt. It can be readily 
shown that 



In conclusion, we have developed an exact and com- 
pletely general reformulation of non-relativistic quantum 
dynamics for pairwise interacting particles by decompos- 
ing the density operator of the A'^-body problem as an av- 
erage of tensor products of stochastic 1-body density op- 
erators. We showed that during the stochastic evolution 
the norm, Hermiticity, and the positivity of the stochastic 
1-body density operators are exactly preserved. Finally 
we explained how the wave vector and hence all quan- 
tum mechanical quantities of interest can be recovered. 
A longer manuscript which presents details of the deriva- 
tion and proves uniqueness of the decomposition is in 
preparation^^. Our equations may have applications as 
a computational technique for chemistry and solid state 
problems, and might also prove useful in themselves as a 
theoretical tool for predicting quantities such as fluctua- 
tions in single molecule fluorescence [l^. 
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*^(0)|Hjv|$^(i')> " *^(* (0)|*^(i')> 
(*^(0)|$^(t')) 



d<'20) 



where the derivative with respect to t' is best computed 
via fast Fourier transform in practice. Since we have a 
closed formula for the phase the exact wave function for 
distinguishable particles reads 



|vI.^(t)) = |$^(i))expH0(i)] 



(21) 



Note that the result (|21|l is itself a sum of tensor products. 
From this expression one may then impose appropriate 
symmetries for indistinguishable particles by operating 
with symmetrization and anti-symmetrization operators. 
Thus the wave vector can be recovered from the density 
operator and hence it is possible to compute all quantum 
mechanical quantities using our approach. It is worth 
noting that average quantities such a,s A — Tr{Ap^} are 
easily computed with our equations but cannot be easily 
obtained with other approaches |^ 0,0. Eigenspectra are 
readily calculated by the Fourier transform of the auto- 
correlation function (*^(0)|*^(t)) using the formula 

I{E) = {^^{0)\d{E-nN)\^''m 
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